# Required fxns: `strs.mean` and `strs.total`

towns=read.csv('http://www.webpages.uidaho.edu/%7Erenaes/Data/tvtowns1.csv',sep=',',header=T)
boxplot(towns$Hours.TV~towns$Town,main='Hours of TV by Town')
Ni=c(93,155,62); ni=c(12,20,8)

strsmean=strs.mean('towns',towns$Hours.TV,towns$Town,ni,Ni)

strstot=strs.total('towns',towns$Hours.TV,towns$Town,ni,Ni)

# use individual vhats from strs.mean and calculate:
# ybar1-ybar2 +/- 2*sqrt(vhat1/n1+vhat2/n2)

# compare TownA and TownB 
# if 0 is in the interval, then they are not significantly different
ybar1=strsmean$ybari[2]
ybar2=strsmean$ybari[3]
vhat1=strsmean$s2i[2]; n1=ni[2]
vhat2=strsmean$s2i[3]; n2=ni[3]
est=ybar1-ybar2; est
bound=2*sqrt(vhat1/n1+vhat2/n2); bound